load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  cs = "dt150" 

  load "loadit" 

  plotname = "figure_MIC_MTB_" + cs 
  plotform = "eps"

  fna = "pc_LP_" + cs + ".nc"
  fnb = "pc_LP_" + cs + ".nc"
  fnc = "pc_LP_" + cs + ".nc"

  fla = addfile(fna,"r")
  flb = addfile(fnb,"r")
  flc = addfile(fnc,"r")

  xa = fla->NITEND_MIC
  xb = flb->NITEND_MTB
  xc = flc->NITEND_MAC

  nt   = dimsizes(xa(:,0,0,0))
  nlev = dimsizes(xb(0,:,0,0))

  v1 = xa(:,:,0,0)
  v2 = xb(:,:,0,0)
  v3 = xc(:,:,0,0)

  vp = new((/3,nlev,nt/),"float")


 wks = gsn_open_wks(plotform,plotname) 
 gsn_define_colormap(wks,"BlueYellowRed") ; choose colormap

; Individual plot resources
  
  res = True
  res@gsnDraw                = False
  res@gsnFrame               = False

  res@cnFillMode = "RasterFill"
 
  res@lbBoxMinorExtentF     = 0.15    
  res@lbOrientation         = "Vertical"
  res@lbLabelsOn            = True
  res@lbLabelFontHeightF    = 0.02
  ;;res@lbLabelAutoStride     = True
  res@lbLabelStride  = 1
  res@lbTitleFontHeightF    = 0.015
  res@tiXAxisFontHeightF    = 0.02
  res@tiXAxisString         = "Days" 
  res@tiYAxisFontHeightF    = 0.02
  res@tiYAxisString         = "Pressure (mb)"
   
;  res@tmXBMinorOn           = False
  res@tmYLMinorOn           = False
  res@tmXTMinorOn           = False
  res@tmYRMinorOn           = False
  res@tmXTOn                = False
  
  ;;res@gsnRightString        = "" 
  ;;res@gsnLeftString         = "" 
  res@gsnStringFontHeightF  = 0.02
  res@gsnCenterStringFontHeightF  = 0.03
  
  res@cnFillOn             = True          ; turn on color
  res@gsnSpreadColors      = True     ; use full colormap
  res@cnLinesOn            = False    ; no contour lines
 
  res@gsnSpreadColors      = True          ; use full range of colormap
  
  ;;if (cmin.ne.cmax) then 
  ;;  res@cnLevelSelectionMode = "ManualLevels"
  ;;  res@cnLevelSpacingF      = cgap
  ;;  res@cnMinLevelValF       = cmin
  ;;  res@cnMaxLevelValF       = cmax
  ;;end if
  res@vpWidthF             = 0.8          ; change aspect ratio of plot
  res@vpHeightF            = 0.3
  ;;res@trYReverse           = True
  
  v2!0 = "time"  ;Reconstruct meta data as the time dimension has changed.
  v2!1 = "lev"
  
  v2&time = xa&time
  v2&lev  = xa&lev
  v2&lev@units = "hPa"
  v2&lev@long_name = "Pressure Levels"

  
  v1!0 = "time"
  v1!1 = "lev"
  
  v1&time = xa&time
  v1&lev  = xa&lev
  v1&lev@units = "hPa"
  v1&lev@long_name = "Pressure Levels"

  v3!0 = "time"
  v3!1 = "lev"

  v3&time = xa&time
  v3&lev  = xa&lev
  v3&lev@units = "hPa"
  v3&lev@long_name = "Pressure Levels"

  res@tmYRMode             = "Automatic"          ; turn off special labels on right axis

  plot = new(4,graphic)   


  res@cnLevelSelectionMode = "ExplicitLevels"    			; set explicit contour levels
  res@cnLevels    = (/ -100.,-50.,-20.,-10.,-5.,-2.,2.,5.,10.,20.,50.,100./) 

  res@gsnLeftString = "a) MIC" 
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 

  plot(0) = gsn_csm_pres_hgt(wks,v1({lev|100:1000},time|:),res)

  res@cnLevels    = (/ -100.,-50.,-20.,-10.,-5.,-2.,2.,5.,10.,20.,50.,100./) 
  res@gsnLeftString = "b) MTB" 
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 
  plot(1) = gsn_csm_pres_hgt(wks,v2({lev|100:1000},time|:),res)

  delete(res@cnLevels) 
  ;;v3 = v1 
  ;;v3 = v2 - v1 
  ;;res@cnLevels    = (/ -5.e6,-2.e6,-1.e6,1.e6,2.e6,5.e6 /) / 1.e9 

  res@gsnLeftString = "c) MAC" 
  res@gsnRightString = "# L~S~-1~N~ s~S~-1~N~" 
  res@gsnRightString = "#/kg/s" ;; L~S~-1~N~ s~S~-1~N~" 
  res@cnLevels    = (/ -100.,-50.,-20.,-10.,-5.,-2.,2.,5.,10.,20.,50.,100./) 
  plot(2) = gsn_csm_pres_hgt(wks,v3({lev|100:1000},time|:),res)


  v4 = v1 
  v4 = (v2 - v1)/(v1+v2+1.e-20) * 100. 
  ;;print(v2 + "    " + v1)

  ;;v4 = where((v1+v2).gt.1.,v4,-999) 
  ;;v4@_FillValue = -999 

  delete(res@cnLevels) 

;;  delete(res@cnLevelSelectionMode) 
  res@cnLevels = (/ -200., -100.,-50.,-20.,-10.,-5., 0., 5.,10.,20.,50.,100.,200./) 

print("P3") 

  res@gsnLeftString = "Rel. Diff" 
  res@gsnRightString = "%" 
  plot(3) = gsn_csm_pres_hgt(wks,v4({lev|100:1000},time|:),res)


;; draw panel with white space added
 resP                 = True
 resP@gsnPanelYWhiteSpacePercent = 5
 resP@gsnPanelXWhiteSpacePercent = 5
 gsn_panel(wks,plot,(/3,1/),resP)


end







